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Abstract 

Many modern statistical applications involve inference for complex stochastic models, 
where it is easy to simulate from the models, but impossible to calculate likelihoods. 
Approximate Bayesian computation (ABC) is a method of inference for such models. 
It replaces calculation of the likelihood by a step which involves simulating artificial 
data for different parameter values, and comparing summary statistics of the simu- 
lated data to summary statistics of the observed data. Here we show how to construct 
appropriate summary statistics for ABC in a semi-automatic manner. We aim for 
summary statistics which will enable inference about certain parameters of interest to 
be as accurate as possible. Theoretical results show that optimal summary statistics 
are the posterior means of the parameters. While these cannot be calculated analyti- 
cally, we use an extra stage of simulation to estimate how the posterior means vary as 
a function of the data; and then use these estimates of our summary statistics within 
ABC. Empirical results show that our approach is a robust method for choosing sum- 
mary statistics, that can result in substantially more accurate ABC analyses than the 
ad-hoc choices of summary statistics proposed in the literature. We also demonstrate 
advantages over two alternative methods of simulation-based inference. 

Keywords: Indirect Inference, Likelihood-free inference, Markov chain Monte Carlo, 
Simulation, Stochastic Kinetic Networks 
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1 Introduction 



1.1 B ackgr ound 



Many modem statistical applications involve inference for stochastic models given 
partial observations. Often it is easy to simulate from the models but calculating 
the likelihood of the data, even using computationally-intensive methods, is 
impracticable. In these natural approach to inference is to use simulations 

from the model for different parameter values, and to compare the simulated data 
sets with the observed data. Loosely, the idea is to estimate the likelihood of a given 
parameter value from the proportion of data sets, simulated using that parameter 
value, that are 'simi l ar to' the observed data. This idea dates back at least as far as 
Diggle and GrattonI (11984 ■ 

Note that if we replace 'similar to' with 'the same as' (see e.g. 



Tavare et al 



19971 ) 



then this approach would give an unbiased estimate of the likelihood; and 
asymptotically as we increase the amount of simulation we get a consistent estimate. 
However in most applications the probability of an exact match of the simulated 
data with the observed data is negligible, or zero, so we cannot consider such exact 
matches. The focus of this paper is how to define 'similar to' for these cases. 

In this paper we focus on a particular approach: approximate Bayesian computation 
(ABC). This approach combines an estimate of the likelihood with a prior to 
produce an approximate posterior, which we will refer to as the ABC posterior. The 
use of ABC initially became popular within population genetics, where simulation 
from a range of p opulation genetic models is possible using the coalescent 
(iKingmanl . Il982[ l. but where calculating likel ihoods is imprac t icable for realistic 



Pritchard et al. 



fll999h who looked at 



sized data sets. The first use of ABC was by 
inference about human demographic history. Fur t her a pplications include inference 
for recombination ra tes (jPadhukasahasram et al.l . 120061 ). evolution of pathogen s 



(jWilson et al.l . l2009l ) and evolution of protein networks ( iRatmann et al. 



20091). Its 



increasing importance can be seen by current r ange of application of ABC, which 



has r ecently been applied wit hin epidemiologv (jMcKinley et al. 



2006 ). inference for extremes flBortot et al 



2009 



20071) ■ dy namical systems ( iToni et al. 



anaka et al. 



2009) and Gibbs random fields f Grelaud et al. . 20091 ) amongst many others. Part of 
the appeal of ABC is its flexibility, it can e asily be applied t o any model for which 
forward simulation is possible. For example IWegmann et al.l ( l2009l ) state that ABC 
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'should allow evolutionary geneticists to reasonably estimate the parameters they 
are really interested in, rather than require them to shift their interest to problems 
for which full-likelihood solutions are available'. Recently there has been software 



developed to help implement ABC withi n population genetics ( ICornuet et al. 



Lopes et al.l . l2009l ) , and systems biology (ILiepe et al. 



2008 



20101 



1.2 ABC Algorithms and Approximations 

Consider analysing n-dimensional data Yobs' have a model for the data, which 
depends on an unknown p-dimensional parameter, 9. Denote the probability density 
of the data given a specific parameter value by 7r(y|^), and denote our prior by 7r(0). 
We assume it is simple to simulate Y from vr(y|6') for any 6*, but that we do not 
have an analytic form for 7r(y|6'). 

We define the ABC posterior in terms of (i) a function S{-) which maps the 
n-dimensional data onto a d-dimensional summary statistic; (ii) a density kernel 
K{x) for a d- dimensional vector x, which integrates to 1; and (iii) a bandwidth 
h > 0. Let s^j^g = 'S'(yQj^g). If we now define an approximation to the likelihood as 

^^(^1^ = / Hy\0)Ki{Siy)-s^^J/h)dy, 
then the ABC posterior can be defined as 

7rABc(^|Sobs) 0^ ^(^)p(^|Sobs)- (1) 

The idea of ABC is that the ABC posterior will approximate, in some way, the true 
posterior for 6 and can be used for inference about 6. The form of approximation 
will depend on the choice of S{-), K{-) and h. For example, if S{-) is the identity 
function then we can view p(^|sQ]gg) as a kernel density approximation to the 
likelihood. If also /i — ?■ 0, then the ABC posterior will converge to the true posterior. 
For other choices of S{-), the kernel is measuring the closeness of y to yQj^g just via 
the closeness of 5'(y) to s^^^g. The reason for considering ABC is that we can 
construct Monte Carlo algorithms which approximate the ABC posterior and only 
require the ability to be able to simulate from 7r(y|6'). 

For simplicity of the future exposition, we will assume that maxK{x.) = 1. This 

assumption imposes no restriction, as if Kq{x.) is a density kernel, then so is 

i^(x) = h^'^KQ^x/ho) for any ho > 0. Thus we can choose Hq so that ma,xK{x) = 1. 
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Note that a bandwidth A for kernel -R'o(x) is equivalent to a bandwidth h = Xh^ for 
K{x.), so the value of ho just redefines the units of the bandwidth h. 

One algorithm for approximating the ABC posterior ([1]), based upon importance 
sampling, is given by Algorithm [TJ (For alternati ve importance sampling 
approaches, based on seq uential Monte Carlo, see 



Beaumont et al. 



2009 



Sisson et al. 



2007, 



20091 ). A standard importance sampler for approximating the 



ABC posterior would repeatedly simulate a paramater value, 6', from a proposal 
g{6) and then assign it an importance sampling weight 



9(0') 



niy\e')K{{S{y) - s^^^}/h)dy 



(2) 



This is not possible here, as vr(y|^) is intractable, so Algorithm [T] introduces an 
extra Monte Carlo step. This step first simulates Ygj^^ from 7r{y\6') and then 
accepts this value with probability -^'^({'S'(ygjj^) — s^j^gj/Zi). Accepted values are 
assigned weights n{6') / g{6'). The key thing is that the expected value of such 
weights is just ([2]), which is all that is required for this to be a valid importance 
sampling algorithm targeting ([1]). 

The output of Algorithm [1] is a weighted sample of 9 values, which approximates the 
ABC posterior of 9. Many of the weights will be 0, and in practice we would remove 
the corresponding 9 values from the sample. A specific case of Algorithm [1] occurs 
when g{9) = it (9), and then we have a rejection sampling algorithm. The most 
common implementation of ABC has a deterministic accept-reject decision in step 
3, and this corresponds to K{-) being the density of a uniform random variable - 
the support of the uniform random variable defining the values of s — s^^^g which 
are accepted. 

An alternative Monte Carlo procedure for implementing ABC is based on MCMC 



(iMarjoram et al.l . 120031 : iBortot et al.l . 120071 ) , and given in Algorithm |2l Both this 
and Algorithr n [J target th e same ABC posterior (for a proof of the validity of this 
algorithm see 



Sisson et al. 



2010[ ). 



Good implementation of ABC requires a trade-off between the approximation error 
between the ABC posterior and the true posterior, and the Monte Carlo 
approximation error of the ABC posterior. The last of these will also be affected by 
the algorithm used to implement ABC, and the specific implementation of that 
algorithm. We will consider both of these approximations in turn. 

To understand the former approximation, consider the following alternative 
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Importance Sampling ABC 
Input: A set of data yQ}-,g, and a function S{-). 

A density kernel A"(-), with maxi^(x) = 1 and a bandwidth h> 0. 
A proposal density g{9); with g{9) > if 7r{6) > 0. 
An integer TV > 0. 

Initialise: Define s^^^ — 'S'(yQ^g). 

Iterate: For i ^ 1, . . . , N 

1. Simulate di from g{0). 

2. Simulate ygjj^^ from 7r(y|0,;), and calculate s = '^'(ysim)- 

3. With probability K{{s — SQ^^)/h set Wi = TT{9i)/ g{9i), otherwise set Wi = 0. 
Output: A set of parameter values {^?i}fli and corresponding weights {wilflj^. 



Algorithm 1: Importance (and rejection) sampling implementation of ABC. 

expression of the ABC posterior as a continuous mixture of posteriors. Consider the 
random variable S = «S'(Y), which is just the summary statistic of the data. Denote 
the posterior for 6 given S = s by 7r(6'|s), and the marginal density for S by 
7r(s) := /7r(s|^)7r(^)d^. Then, 

...c(»|s„b,) = / ^(s).(»|s)ds. Where ^(s) = j^'^Z^'/^^^ ^ (3) 

The mixing weight for a given value of S is just the conditional density for such a 
value of S given acceptance in step 3 of Algorithm [TJ If /i is small then 
'^abc{&\^q})s) ~ ^(^l^obs)' Thus we can split the approximation of the posterior by 
ABC into the approximation of 7i"(^|yobs) '''"(^l^obs)' further 
approximation of -7r(0|sQ|^g 

) by ^ABc(^|sQ]3g). The former is controlled by the choice 
of S{-), the latter by K{-) and h. 

Now consider the Monte Carlo approximation within importance sampling ABC. 
Consider a scalar function h{6), and define the ABC posterior mean of h{9) as 

EABc(/i(^)|Sobs)= / /^WvrABc(^|Sobs)d^- 
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MCMC ABC 



Input: A set of data yQj^g, and a function S{-). 

A density kernel K{-)^ with niaxi^(x) = 1 and a bandwidth h> 0. 
A transition kernel g{-\-)- 
An integer N > 0. 

Initialise: Define s^^^ = '^'(yobs)' ^^"^ choose or simulate Oq and Sq. 

Iterate: For i = 1, . . . , N 

1. Simulate 6 from g{9\9i^i). 

2. Simulate Ygjjj^ from 7r(y|0), and calculate s = 5'(ygjjj^)- 

3. With probability 

. r A-((s - s^bs)/^) 7r(%(g,,i|g) ] 
miri s. I > 

1 '/v((s,_i-s„bs)//*)4&,_i)5(^|e.-i)J 
accept 6* and set 9i ~ 9 and S; = s; otherwise set 9i = 9i-i and = Si_i. 

Output: A set of parameter values {9i}f^i. 



Algorithm 2: MCMC sampling implementation of ABC. 
Assuming this exists then, by the law of large numbers, as — )■ oo 

%^^^E.ecW)|Sobs), (4) 

where convergence is in probability. For rejection sampling, for large A^ the variance 
of this estimator is just 

VarABc(/^(^)|Sobs) 
^^acc 

where the numerator is the ABC posterior variance of h{9)^ and the denominator is 
A'acc = A^ / p(^|Sobs)'''"(^)'^^' expected number of acceptances. 
Similar calculations for both importance sampling in general and for the MCMC 
version of ABC are given in Appendix A. In all cases the Monte Carlo error depends 
on J p(6'|sQ]gg)7r(6')d6', the average acceptance probability of the rejection sampling 
algorithm. The following Lemma characterises this acceptance probability for small 
h in this case of continuous summary statistics (similar results can be shown for 
discrete summary statistics). 
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Lemma 1 Assume that either (i) the marginal distribution for the summary 
statistics, 7r(s) is continuous at s = s^^^ and that the kernel A'(-) has finite support; 
or (a) 7r{s) is continuously diff'erentiable, \d7T{s)/dsi\ is bounded above for all i, and 
J |xj|ii'(x)(ix is bounded for all i. Then, in the limit as h ^ 

J pie\s^f^^Me)dB = 7r(s^g/^'^ + (/i^) , (6) 

where d is the dimension ofs^^^. 

Proof: See Appendix B. 

This result gives insight into how S{-) and h affect the Monte Carlo error. To 
minimise Monte Carlo error, we need h'^ to not be too small. Thus ideally we want 
S{-) to be a low-dimensional summary of the data that is sufficiently informative 
about 9 so that 7r(6'|sQ|^g) is close, in some sense, to 7r(6'|yQ|^g). The choice of h 
affects the accuracy of ttabc in approximating ^T{9\s^^^^, but also the average 
acceptance probability ([6]), and hence the Monte Carlo error. 



1.3 Approach and Outline 

Previous justification for ABC has, at least informally, been based around ttabc 
approximating T^iG\yQ\^^) globally. This is possible if Yq]^^ is low-dimensional, and 
S{-) is the identity function, or if we have low-dimensional sufficient statistics for 
the data. In these cases we can control the error in the ABC approximation by 
choosing h sufficiently small. In general applications this is not the case. Arguments 
have been made about tryi ng to choose approximate sufficient statistics (see e.g. 



Joyce and MarjoramL l2008l ). However the definition of such statistics is not clear, 
and more importantly it is difficult or impossible to construct a general method for 
finding such statistics. 

We take a different approach, and weaken the requirement for vtabc to be a good 
approximation to niOlyQ-^^. We argue for vtabc to be a good approximation solely 
in terms of the accuracy of certain estimates of the parameters. We also would like 
to know how to interpret the ABC posterior and probability statements derived 
from it. To do this we consider a property we call calibration, which we define 
formally below. If ttabc is calibrated, then this means that probability statements 
derived from it are appropriate, and in particular that we can use ttabc to quantify 
uncertainty in estimates. We argue that using such criteria we can construct ABC 
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posteriors which have both good inferential properties, and which can be estimated 
well using Monte Carlo methods such as in Algorithm [TJ 

In Section [2] we define formally what we mean be calibration and accuracy. 
Standard ABC is not calibrated, but a simple modification, which we call noisy 
ABC, is. We further show that if the bandwidth, h, is small then standard ABC is 
approximately calibrated. Theoretical results are used to produce recommendations 
about when standard ABC and when noisy ABC should be used. Then, for a 
certain definition of accuracy we show that the optimal choice of summary statistics 
S{Y) are the true posterior means of the parameters. Whilst these are unknown, 
simulation can be used to estimate the summary statistics. Our approach to doing 
this is described in Section |3l together with results that show the advantage of ABC 
over directly using the summary statistics to estimate parameters. Section H] gives 
examples comparing our implementation of ABC with previous methods. The paper 
ends with a discussion. 

2 Calibration and Accuracy of ABC 

Firstly we define what we mean by calibration, and introduce a version of ABC that 
is calibrated. We show how the idea of calibration can be particularly important 
when analysing multiple data sets. We then discuss our definition of accuracy of the 
ABC posterior, and how this can be used to guide the choice of summary statistic 
we use. 



2.1 Calibration and Noisy ABC 



Consider a subset of the parameter space A.. For given data Ygj^g, the ABC 
posterior will assign a probability to the event 9 ^ A 




For a given probability, q, consider the event Eg{A) that PrABc(^ ^ -^l^obs) ~ ^' 
Then the ABC posterior is calibrated if 



Pt{9 e A\E,{A)) = Q- 



(7) 
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The probability is then defined in terms of the density on parameters and data 
given by our prior and hkehhood model, 7r(^)7r(y|6'), but ignoring any Monte Carlo 
randomness. Statement states that under repeated sampling from the prior, 
data, and summary statistics, events assigned probability q by the ABC posterior 
will occur with probability q. A consequence of calibration is that the ABC 
posterior will appropriately represent uncertainty in the parameters: for example we 
can construct appropriate credible intervals. That is calibration means we can use 
the ABC posterior as we would any standard posterior distribution. 

Standard ABC posteriors ([3]) are not calibrated in general. Instead we introduce the 
idea of noisy ABC which is calibrated. Noisy ABC involves defining summary 
statistics which are random. A noisy ABC importance sampling algorithm is 
obtained by changing the initialisation within Algorithm [1] or Algorithm |2l details 
are given in Algorithm [31 The resulting ABC posterior for a given y^]^^ is random, 
and the definition of the probability in ([7]) needs to account for this extra 
randomness. 

Noisy ABC 

Replace the initialisation step in Algorithm [1] and [2] with: 
Initialise: Simulate x from A'(x). Define s^^g = "^(yobs) + 



Algorithm 3: Change of Algorithm [T] or [2] for noisy ABC. 

Theorem 1 Algorithmic produces an ABC posterior that is calibrated. 

Proof: The ABC posterior derived by Algorithm [3] is vrABc(^|sQi3g) where s^j^g is 
related to the data by YqI^^ by 

Sobs = *^(yobs) + ^^' (8) 
and X is the realisation of a random variable with density K(x.). However the 
definition of 7rABc(^|sQ]3g) is just that of the true posterior for 9 given data s^^^g 
generated by ([S]). It immediately follows that this density is calibrated. □ 



A related idea is used by IWilkinsonI (j2008[ l who shows that the ABC posterior is 



equivalent to the true posterior under an assumption of appropriate model error. In 
the limit as /i — )■ noisy ABC is equivalent to standard ABC, we discuss the links 
between the two in more detail in Section |2T 
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2.2 Inference from Multiple Data Sources 

Consider combining data from m independent sources, Yoi^g) • • • 'yg^bs' possible 
to use individual ABC analyses for each data set, and then combine these 
inferences. One sequential approach is to use the ABC posterior after analysing the 
2th data set as a prior for analysing the {i + l)st data set. Algorit hms for such an 
approach are a special case of those discussed in Wilkinson f 201ll ). 



One consequence of calibration is that such inferences will be well behaved in the 
limit as m gets large. To see this define the ABC approximation to the ith data set 
as 

^^^l^obs^' 

and note that the above approach is targeting the follow ABC posterior 

m 

w(«isx.....ib',)«M»)np(''Ni)- 

i=l 

If we use noisy ABC, where s^j^^ is random centered on 'S'(y^*|^g), then by same 
argument as for Thereom [T] this ABC posterior will be calibrated. Furthermore, we 
have the following result 

Theorem 2 Let 6q be the true parameter value. Consider noisy ABC, where 
^obs ~ ^(^obs^ ~^ '^^^''"^ ^ drawn from K{-), then the expected noisy-ABC 
log-likelihood, 

E{\og[p{e\S^f^^)]}= J J\og[p{e\S{y) + h^)My\eo)K{^)dyd^, 
has its maximum at 6 = 6q. 

Proof: Make the change of variable s^^^g = S(y) + /ix, then 

E{log[p(^|Sobs)]} = ;^ / / log[p(e|Sobs)]^(y|^o)i^([Sobs-S(y)]//^)dydSobs- 
Now by definition p{6\s^^-^^ = f 7r(y|6')Ji ([s^^g — S(y)]//i)dy, thus we get 

E{log[p(^|Sobs)]} = \jj log[p(^|Sobs]P(^l^obs)d^obs- 

and by Jensen's inequality this has its maximum at ^ = 9q. □ 
he importance of th i s resu lt, is that under the standard regularity conditions 



(IBernardo and Smithl . Il994l ). the ABC posterior will converge onto a point mass on 
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the true parameter value as m — )■ oo. By comparison, if we use standard ABC, then 
we have no such guarantee. As a simple example assume F*-*-* are independent 
identically distributed from a Normal distribution with mean and variance o"^, and 
our kernel is chosen to be Normal with variance r < cr^. If we use standard ABC, 
the ABC posterior given y'^^\ . . . , y^"^^ will converge to a point mass on cr^ — r as 
m oo. 

We look empirically at this issue in Section 14.31 

2.3 Accuracy and Choice of Summary Statistics 

Calibration itself is not sufficient to define a sensible ABC posterior. For example 
the prior distribution is always calibrated, but will not give accurate estimates of 
parameters. Thus we also want to maximise accuracy of estimates based on the 
ABC posterior. We will define accuracy in terms of a loss function for estimating 
the parameters. A natural choice of loss function is quadratic loss. Let 6q be the 
true parameter values, and 6 an estimate. Then we will consider the class of loss 
functions, defined in terms of a p x p positive-definite matrix A, 

Li9o,9;A) = i9o-9fA{9o-9). (9) 

We now consider implementing ABC to minimise this quadratic error loss of 
estimating the parameters. We consider the limit oi h ^ 0. This will give results 
that define the optimal choice of summary statistics, S{-). 

For any choice of weight matrix A that is of full-rank, the following theorem shows 
that the optimal choice of summary statistics is 5'(yQ]^g) = E{9\-y^^^), the true 
posterior mean. 

Theorem 3 Consider a p x p positive- definite matrix A of full rank. Given 
observation y^^g, let S be the true posterior variance for 9. Then 

(i) The minimal possible quadratic error loss, E{L{9^9] A)\y occurs when 
9 = E{9\'ygjjg) and is trace{AIl). 

(a) If Siygjjg) = E{9\'y^^g) then in the limit as h ^ the minimum loss, based on 
inference using the ABC posterior, is achieved by 9 = E^^gci^l^ obs^ ■ 
resulting expected loss is trace{AT,). 
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Proo f: Part (i) is a standard result of Bayesian statistics (iBernardo and Smith . 



19941 ) ■ For (ii) we just need to show that in the hmit as — ?■ 0, 

EABc(^|sQigg) = E(^|yQ]^g). By definition in the hmit h s^j^^g = 5'(yQ]gg) with 

probabihty 1, and vrABc(^|sQi3g) = 7r(6'|sQ]^g). Furthermore 

EABc(^lsobs) = / 07T{e\s^^^)de, 

= 11 ^7r(^|y)7r(y|s^bs)dyd^' 

where 7r(y|sQj^g) is the conditional distribution of the data y given the summary 

statistic s^Y)s- Finally by definition, all y that are consistent with s^j^^ satisfy 

J ^7r(^|y)d^ = E(^|yQ]gg), and hence the result follows. □ 

Use of square error loss leads to ABC approximations that attempt to have the 
same posterior mean as the true posterior. Using alternative loss functions would 
mean matching other features of the posterior: for example absolute error loss 
would result in matching the posterior medians. It is possible to choose other 
summary statistics that also achieve the minimum expected loss. However any such 
statistic with dimension d > p will cause larger Monte Carlo error (see Lemma [1] 
and the discussion below). 



2.4 Comparison of standard ABC and noisy ABC 

Standard and noisy ABC are equivalent in the limit as /i — > 0, with the ABC 
posteriors converging to E(^|S(yQ]2g)). We can further quantify the accuracy of 
estimates based on standard or noisy ABC for h ^ 0. For noisy ABC we have the 
following result 

Theorem 4 Assume condition (i) of LemmaUi that 7r(i?(^|y^^g)) > 0, and the 
kernel K{-) corresponds to a random variable with mean zero. If 
S(yQ^g) = E{9\y ^^^^ then for small h the expected quadratic loss associated with 
9 = EABc{d\sQi)g) is 

E{L{e, 9- A)\y^Q = traceiAE) + J x^AxK{x)dx + o{h^). 

Proof: The idea is that EABcl^lsQj^g) = s^j^g + o(/i); and the squared error loss 
based on ^ = s^j^g is just trace(AS) + h"^ J x^A:x.K{x.)d'x. See Appendix C. 
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A similar result exists for standard ABC (jPrangld . |2010[ ). which shows that in this 



case 

,4 



E{L{e, 9; A) ly^^g) = trace(AS) + 0{h^ 

Extensions of these results can be used to give guidance on the choice of kernel. For 
noisy ABC they suggest a uniform kernel on the ellipse x'^Ax < c, for some c, for 
standard ABC they also suggest a uniform kernel on an ellipse, but the form of the 
ellipse is difficult to calculate in practice. We do not give these results in more 
detail, as in practice we have found the choice of kernel has relatively little effect on 
the accuracy of either ABC algorithm. 

These two results also give an insight into the overall accuracy of a Monte Carlo 
ABC algorithm. For simplicity consider a rejection sampling algorithm. The Monte 
Carlo variance is inversely proportional to the acceptance probability. Thus using 
Lemma [1] we have that the Monte Carlo variance is 0{N~^h~'^), where is the 
number of proposals. The expected quadratic loss based on estimates from the ABC 
importance sampler will been increased by this amount. Thus for noisy ABC we 
want to choose h to minimise 

Co 



trace(AS) + J x^y4xfr(x)dx 



for some constant Cq. This gives that we want h = 0{N~^^^'^~^'^^), and the overall 
expected loss above trace(74S) would then decay as A^~2/(2+d)_ standard ABC a 
similar argument give h = 0{N~^^^^'^^^), and the overall expected loss would then 
decay as Ar-4/(4+d)_ 

Thus we can see that the choice between using standard ABC and noisy ABC is one 
of a trade-off between accuracy and calibration. Noisy ABC is calibrated, but for 
small h will give less accurate estimates. As such for the analysis of a single data set 
where the number of summary statistics is not too large, and hence h is small, we 
would recommend the use of standard ABC. If we wish to combine inferences from 



ABC analyses of multiple data sets, then in light of the discussion in Section \2.2\ we 
would recommend noisy ABC. For all the examples in Section H] we found that this 
approach worked well in practice. 

Note that possibly the best approach is to use noisy ABC, but to use 
Rao-Blackwellisation ideas to average out the noise that is added to the summary 
statistics. Such an approach would have the guarantee that the resulting expected 
quadratic loss for estimating any function of the parameters would be smaller than 
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that from noisy ABC. However implementing such a Rao-Blackwelhsation scheme 
efficiently appears non-trivial, with the only simple approach being to run noisy 
ABC independently on the same data set, and then to average the estimates across 
each of these runs. 

3 Semi-automatic ABC 

The above theory suggests that we wish to choose summary statistics that are equal 
to posterior means. Whilst we cannot use this result directly, as we cannot calculate 
the posterior means, we can use simulation to estimate appropriate summary 
statistics. 

Our approach is to: 

(i) use a pilot run of ABC to determine a region of non-negligible posterior mass; 

(ii) simulate sets of parameter values and data; 

(iii) use the simulated sets of parameter values and data to estimate the summary 
statistics; and 

(iv) run ABC with this choice of summary statistics. 

Step (i) of this algorithm is optional. Its aim is to help define an appropriate 
training region of parameter space from which we should simulate parameter values. 
In applications where our prior distributions are relatively informative, this step 
should be avoided as we can simulate parameter values from the prior in step (ii). 
However it is important if we have uninformative priors, particularly if they are 
improper. 

If we implement (i), we assume we have arbitrarily chosen some summary statistics 
to use within ABC. In our implementation below we choose our training region as a 
hypercube, with the range for each parameter being the range of that parameter 
observed within our pilot run. Then in (ii) we simulate parameter values from the 
prior truncated to this training region, and for each choice of parameter value we 
simulate an artificial data set. We repeat this M times, so that we have M sets of 
parameter values, each with a corresponding simulated data set. 
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There are various approaches that we can take for (iii). In practice we found that 
using hnear regression, with appropriate functions of the d ata as predictors, is both 



simple and worked well. We al so considered using LASSO (iHastie et al.l . 120011 ) . and 
canonical correlation analysis ( iMardia et al.l . 119791 ) but in general neither of these 
performed better than linear regression (though the LASSO maybe appropriate if 
we wish to use a large number of explanatory variables within the linear model). 

Our linear regression approach involved considering each parameter in term. First 
we introduce a vector- valued function /(■), so that /(y) is a vector of, possibly 
non-linear, transformations of the data. The simplest choice is /(y) = y, but in 
practice including other or different transformation as well may be beneficial. For 
example, in one application below we found /(y) = (y, y^, y'^, y^), that is a vector of 
length 4n that consists of the data plus all second to fourth powers of individual 
data points, produced a better set of summary statistics. 

For the ith summary statistic the simulated values of the ith parameter, 

6^^\ . . . , 6l^^\ are used as the responses; and the transformations of the simulated 

data, f{y^^^), . . . , /(y*-*^^), are used as the explanatory variables. We then fit the 

model 

^, = E(^,|y) + e, = /3j^) + /3«/(y) + e^, 

where ej is some mean-zero noise, using least squares. The fitted function 

+ /3^*'*/(y) is then an estimate of E(6'j|y). The constant terms can be neglected 
in practice as ABC only uses the difference in summary statistics. Thus the ith 
summary statistic for ABC is just /3*-*''/(y). 

Note that our approach of using a training region means that our model for the 
posterior means are based only on parameter values simulated within this region. 
We therefore suggest adapting the ABC run in step (iv) so that the prior is 
truncated to lie within th is training region (a similar idea is used in 
Blum and Francoisl . 2010). This can be viewed as using, weakly, the information we 



have from the pilot ABC run within th e final ABC run, and has links with 



composite likelihood methods (iLindsayl . 119881 ). More importantly it makes the 
overall algorithm robust to problems where E(/3'^*^/(Y)|6'i) is similar for two 
dissimilar values of 6i, one inside the training region and one outside. 

In practice below we use roughly a quarter of our total CPU time on (i) and (ii) and 
half on (iv), with step (iii) having negligible CPU cost. Note that we call this 
semi-automatic ABC as the choice of summary statistics is now based on 
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simulation, but that there are still choices by the user in terms of fitting the linear 
model in step (iii). This input is in terms of the choice of /(y) to use. Note that 
(iii) is now a familiar statistical problem, and that standard model checks can be 
used to decide whether that choice of /(y) is appropriate, and if not, how it could 
be improved. Also note that repeating step (iii) with a different choice of /(y) can 
be done without any further simulation of data, and thus is quick in terms of the 
CPU cost. Furthermore standard model comparison procedures (e.g. using BIG) 
can be used to choose between summary statistics obtained from linear-regressions 
using different explanatory variables. 

A natural question is whether this approach is better than the current approach to 
ABC, where summary statistics are chosen arbitrarily. In our implementation we 
still need to choose summary statistics for step (i), and we also need to choose the 
set of explanatory variables for the linear model. However we believe our approach 
is more robust to these choices than standard ABC. Firstly the choice of summary 
statistics in step (i) is purely to make step (ii) more efficient, and as such the final 
results depend little on this choice. Secondly, when we choose the explanatory 
variables we are able to choose many such variables (of the order of 100s). As such 
we are much more likely to include amongst these some variables which are 
informative about the parameters of interest than standard ABC where generally a 
few summary statistics are used. If many summary statistics are used in ABC, then 
this will require a large value of /i, and will often be inefficient due to the 
accept/reject decision being based not only on the informative summary statistics, 
but also those which are less informative. These issues are demonstrated empirically 
in the examples we consider. 

Our approach has s i milar ities to that of 



Beaumont et al. 



( 120021 ) (see also 



Blum and Francoisl . 120101 ) , where the authors use linear regression to correct the 



output from ABC. The key difference is t hat our approach uses li near regression to 



construct the summary statistics, whereas 



Beaumont et al. 



(I2OO2I ) uses linear 



regression to reduce the error bet ween 7rABc(^|sQ]Qg) and 7r(6'|sQ|^g). In particular the 



method of Beaumont et al 



( I2OO2I ) assumes that appropriate low-dimensional 



summary statistics h ave already been c hosen . We look at differences between our 



approach and that of 



Beaumont et al. 



( I2OO2I ) empirically in the examples. 
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3.1 Why use ABC? 



Our approach involves using simulation to find estimates of the posterior mean of 
each parameter. A natural question is why not use these estimates directly? We 
think using ABC has two important advantages over just using these estimates 
directly. The first is that ABC gives you a posterior distribution, and thus you can 
quantify uncertainty in the parameters as well as getting point estimates. 

Moreover we have the following result. 

Theorem 5 Let 9 = i?(6'|5'(y^^g)). Then for any function g, 

E {l{9, ~9- A) |5(y^g) < E {1(9, ^?(5(y^g); A) |5(y^g) . 

Furthermore, asymptotically as h the ABC posterior mean estimate of 9 is 
optimal amongst estimates based on S{ygjjg). 

Proof: The proof of the first part is t he standard argument that t he mean is the 



optimal estimator under quadratic loss 



Bernardo and Smith 



Jl994[ ). 



The second 



part follows because as /i — )■ the ABC posterior mean tends to 9. □ 

Note that this result states that, in the limit as /?,—)■ 0, ABC gives estimates that 
are at least as or more accurate than any other estimators based on the same 
summary statistics. 

Comparison with Indirect Inference 



Indirect inference (IGourieroux and Ronchetti Il993l ) is a method similar to ABC in 



that it uses simulation from a model to produce estimates of the model's 
parameters. The general procedure involves first analysing the data under an 
approximating model, and estimating the parameters, called auxilary parameters, 
for this model. Then data is simulated for a range of parameter values, and for each 
simulated data set we get an estimate of the auxilary parameters. Finally we 
estimate the true parameters based on which parameter values produced estimates 
of the auxilary parameters closest to those estimated from the true data. (In 
practice we simulate multiple data sets for each parameter value, and get an 
estimate of the auxilary parameters based on these multiple data sets.) The link to 
ABC is that the auxilary parameters in indirect inference are equivalent to the 
summary statistics in ABC. Both methods then use (different) simulation 
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approaches to produce estimates of the true parameters from the values of the 
auxilary paramaters (or summary statistics) for the real data. 

For many approximating models, the auxilary parameters depend on a small set of 
summary statistics of t he da ta; these are called auxilary statistics in 



Heggland and Frigessil (120041 ) . In these cases indirect inference is performing 



inference based on these auxilary statistics. The above result shows that in the limit 
as /i — 0, ABC will be more accurate than an indirect inference method whose 
auxilary statistics are the same as the summary statistic used for ABC. We 
investigate this empirically in Section HI 



4 Examples 



The performance of semi-automatic ABC was investigated in a range of examples: 
independent draws from a complex distribution (Section 14. 2p . a stochastic kinetic 
network for biochemical reactions (Section l4.3p . a partially observed M/G/1 queue 
(Section l4.5p . an ecological population size model (Section l4.4p . and a Tuberculosis 
transmission model (Section 14. 6p . Section \A7i\ describes implementation details that 
are common to all examples. Section 14.31 concerns a dataset for which ABC and 
related methods have not previously been used, and highlights the use of noisy ABC 
in the sequential approach of Section 12.21 The other examples have previous 
analyses in the literature, and we show that semi-automatic ABC compares 
favourably against e xisting r netho ds including indirect inference, the synthetic 
likelihood method of 



WoodI (120101 ) and ABC w ith ad-hoc sumrnary sta tistics, with 



Beaumont et al. 



(l2002h . We also 



or without the regression correction method of 
show that direct use of the linear predictors created during semi-automatic ABC 
can be inaccurate (e.g Section l46l) . Outside Section 1473) noisy ABC runs are not 
shown; they are similar to non-noisy semi-automatic ABC but slightly less accurate. 
The practical details of implementing our method are also explored, in particular 
how the choice of explanatory variables /(y) is made. 



4.1 Implementation details 

Apart from the sequential implementation of ABC in 14.31 ^ill ABC analyses were 
performed using Algorithm [2] with a Normal transition kernel. The density kernel 
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was uniform on an ellipsoid x^Ax < c. This is a common choice in the literature 
and close to optimal for runs using semi-automatic ABC summary statistics as 
discussed in Section 12.41 For summary statistics not generated by our method we 
generally we used A = I. In Section 14.41 a different choice was necessary and is 
discussed there. For ABC using summary statistics from our method, recall from 
Section 12.31 that A defines the relative weighting of the parameters in our loss 
function. In Sections 14.21 and 14.61 the parameters are on similar scales so we used 
A = I. Elsewhere, marginal parameter variances were calculated for the output of 
each pilot run and the means of these [sf, s^, . . .) taken. A diagonal A matrix was 
formed with zth diagonal entry s^'^. 

Other tuning details required by Algorithm [2] are the choice of h, the variance 
matrix of the transition kernel and the starting values of the chain. Where possible, 
these were chosen by manual experimentation or based on previous analyses 
(e.g. from pilot runs). Otherwise they were based on a very short ABC rejection 
sampling analysis. Except where noted otherwise, h was tuned to give an 
acceptance rate of roughly 1% as this gave reasonable results in the applications 
considered. (An alternati ve would be to use computational rnethods that try and 
choose h for each run, see 



Bortot et al. 



2007; 



Ratmann et al. 



20071). 



In the following examples our method is compared against an ABC analysis with 
summary statistic based on the existing literature, referred to as the "comparison" 
analysis. To allow a fair comparison, this uses the same number of simulations as 
the entire semi-automatic method and a lower acceptance rate, roughly 0.5%. 

For simulation studies on multiple datasets, the accuracies of the various analyses 
were compared as follows. The point estimate for each dataset was calculated, and 
the quadratic loss ([9]) of each parameter estimate relative the true parameter value 
was calculated. We present the mean quadratic losses of the individual parameters. 
In tables of results we highlight the smaller quadratic losses (all with 10% of the 
smallest values) by italicising. 



4.2 Inference for g and k Distribution 



The g and k distribution is a flexibly shaped distribution used to model 
non-standard data through a small number of parameters (iHaynesl . Il998l ). It is 
defined by its inverse distribution function (ITOl) . below, but has no closed form 
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density. Likelihoods can be evaluated numerically but this i s costly 
(iRayner and MacGillivray . 



2002 



Drovandi and Pettit 



20091 ). ABC methods are 



attractive because simulation is straightforward by the inversion method. Here we 
use the fact we can calculate the maximum likelihood estimate numerically to also 
compare ABC with a full-likelihood analysis. 

The distribution is defined by: 

1 



F-\x; A,B,c,g,k) = A + B [1 + c 



exp[-gz[x) 



[l + z(x?)''z(x) 



(10) 



1 + exp{—gz{x) 

where z{x) is the xth standard Normal quantile, A and B are location and scale 
parameters and g and k are related to skewness and kurtosis. The final parameter c 
is typically fixed as 0.8, and this is assumed throughout, leaving unknown 
parameter s 6 = (A, B, q, k). The only para. meter restrictions are B > and 
k > -1/2 jRavner and MacGillivravl . 120021 ). 

AUingham et al.l (120091 ) used ABC to analyse a simulated dataset oi n = 10^ 
independent draws from the g-and-k distribution with parameters 6o = (3, 1, 2, 0.5). 
A uniform prior on [0, 10] was used and the summary statistics were the full set of 
order statistics. We studied multiple datasets of a similar form as detailed below. 
Our aim is firstly to show how we can implement semi-automatic ABC in a 
situation where there are large numbers of possible explanatory variables (just using 
the order statistics gives 10^ explanatory variables), and to see how the accuracy of 
se mi-automatic ABC com pares with the use of arbitrarily chosen summary statistics 



m 



AUingham et al.l ( 120091 ). We also air a to look at comparing semi-automatic ABC 

(I2OO2I ). and with indirect 



Beaumont et al. 



with the linear regression correction of 
inference. 

Comparison of ABC methods 

The natural choice of explanatory variables for this problem are based upon the 
order statistics, and also powers of the order statistics. Considering up to the fourth 
power seems appropriate as informally the four parameters are linked to location, 
scale, skewness and kurtosis. However fitting the linear model with the resulting 
4 X 10^ exaplanatory variables is impracticable. As a result we considered using a 
subset of m evenly spaced order statistics, together with up to I powers of this 
subset. To choose appropriate values for m and / we fitted linear models with m 
ranging over a grid of values between 60 and 140, and / ranging between 1 and 4, 
and used BIC (averaged across the models for the four parameter values) to choose 
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an appropriate value for m and /. We then used the summary statistics obtained 
from the hnear model with this value of m and I in the final run of ABC. For 
simplicity we did this for the first data set, and kept the same value of m and / for 
analysing all subsequent data sets. 

Note that using subsets of order statistics has computational advantages, as these 
can be generated efficiently by simulating correspor iding standar d uniform order 
statistics using the exponential spacings method of iRipleyl (119871 ) (page 98) and 
performing inversion by substituting these in (fTOjl . The cost is linear in the number 
of order statistic s requ ired. Our pilot ABC run used the summary statistics from 



AUingham et al.l (120091 ). Fitting the different linear models added little to the 



overall computational cost, which is dominated by the simulation of the data sets at 
the different stages of the procedure. 

Our semi-automatic ABC procedure chose m = 100 order statistics and / = 4. The 
accuracy of the resulting parameter estimates, measured by square error loss across 
implementation of semi-automatic ABC on 50 data sets is shown in Table [H For 
comparison we show results of the lAllingham et al.l ( l2009l ) ABC method, 
implemented to have the same overall computational cost. We also show the 
accuracy of estir nates obtained by post- processing the lAllingham et al.l (120091) 
results using the iBeaumont et al.l (120021 ) regression correction. We could only use 
this regression correction on 48 of the 50 data sets, as for the remaining 2, there 
were too few acceptances in the ABC run for the regression correction to be stable 
(on those two data sets the resulting loss after performing the regression correction 
was orders of magnitude greater than the original ABC estimates). 

Despite the lAllingham et al.l ( 120091 ) analysis performing poorly, and hence 
producing a poor pilot region for semi-automatic ABC; semi-automatic ABC 
appears to perform well. It has losse s that are between a factor of 2 and 100 smaller 
than the method of lAllingham et al.l ( 120091 ) with or without the regression 
correction. Using the regression correction does improve accuracy of the estimates 
for three of the four parameters, but to a lesser extent than semi-automatic ABC. 
For comparison, using the predictors from the linear-regression to directly estimate 
the parameters had similar accuracy for all parameters except g, where the 
linear-predictor's average error was greater by about a third. 

To investigate the effect of the pilot run on semi-automatic ABC, and the choice of 
summary statistics on the regression correction, we repeated this analysis by 
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implementing ABC with 100 order statistics. This greatly improved the 
performance of ABC, showing the importance of the choice of summary statistics. 
The improved pilot run also improves semi-automatic ABC but to a much lesser 
extent. In this case semi-automatic ABC has similar accuracy to the comparison 
ABC run with the regression correction. For a further comparison we also 
calculated the MLEs for each data set numerically. The two best ABC runs have 
mean quadratic loss that is almost identical to that of the MLEs. 





A 


B 


g 


k 


Allingham et al 


0.0059 


0.0013 


3.85 


0.00063 


Allingham + reg 


0.00040 


0.0017 


0.28 


0.00051 


Semi-Automatic ABC 


0.00016 


0.00056 


0.044 


0.00023 


Comparison 


0.00025 


0.00063 


0.0061 


0.00041 


Comparison + reg 


0.00016 


0.00055 


0.0014 


0.00015 


Semi-Automatic ABC 


0.00015 


0.00053 


0.0014 


0.00015 


MLE 


0.00016 


0.00055 


0.0013 


0.00014 



Table 1: Mean quadratic losses of various ABC analyses of 50 g-and-k datascts with parameters 
(3,1,2,0.5). The first three rows are based on using the Allingham et al. summary statistics in 
ABC, and in the ABC pilot run for semi-automati c ABC. The nex t three rows use just 100 evenly 
spaced order statistics. Results based on using the lBeaumont et al.l (|2002l ) regression correction are 
denoted reg. For Allingham + reg, we give the mean loss for just 48 of the 50 data sets. For the 
remaining 2 data sets the number of ABC acceptances was low sa 200, and the regression correction 
was unstable. For comparison we give the mean quadratic loss of the true MLEs. 



Comparison with Indirect Inference 

Theorem O shows that asymptotically ABC is at least as accurate as other 
estimators based on the same summary statistics. We tested this with a comparison 
against indirect inference (see Section [3TTi) . The semi-automatic ABC analysis was 
repeated, and to give a direct comparison, its summary statistics were used as the 
indirect inference auxiliary statistics. 

Initial analysis showed that which method is more accurate depends on the true 
value of 6 and in particular the parameter g; this is illustrated by Figure [1] 
Therefore we studied datasets produced from variables g values; we drew 50 g values 
from its prior and for each simulated datasets of n g-a.nd-k draws conditional on 
e = (3, 1, g, 0.5) for n = 10^, 10=^ and 10^ 
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Each semi-automatic ABC analysis used a total of 3.1 x 10^ simulated data sets. 
Indirect inference was roughly tuned so that the total number of simulations 
equalled this. (Similar results can be obtained from indirect inference using many 
fewer simulations, and indirect inference is thus a computationally quicker 
algorithm). Mean losses are given in Table IH showing that while the methods 
perform similarly for n = 10^, ABC is more accurate for smaller n. 

More detail is given in Figure [U which plots the true against estimated g values. Of 
particular interest is the n = 100 case where the g parameter is very hard to identify 
for g > 3. It is over this range that ABC out-performs indirect inference most 
clearly, with estimates from indirect inference being substantially more variable 
than those for ABC. 









A 


B 


g 


k 






Pilot 


0.0003 


0.0008 


1.7 


0.0004 


n = 


10^ 


Indirect inference 


0.0003 


0.0022 


0.082 


0.0063 






Semi-automatic ABC 


0.0001 


0.0005 


0.059 


0.0002 






Pilot 


0.0031 


0.014 


4.6 


0.0073 


n = 


103 


Indirect inference 


0.0066 


0.014 


0.83 


0.0053 






Semi-automatic ABC 


0.0012 


0.0094 


0.51 


0.0042 






Pilot 


0.0089 


0.039 


4.8 


0.057 


n = 


102 


Indirect inference 


0.018 


0.059 


5.5 


0.067 






Semi-automatic ABC 


0.0075 


0.046 


3.5 


0.040 



Table 2: Mean quadratic losses of semi-automatic ABC and indirect inference analyses of 50 
g-and-fc datascts with variable g parameters. 



4.3 Inference for Stochastic Kinetic Networks 

Stochastic kinetic networks are used to model biochemical networks. The state of 
the network is determined by the number of each of a discrete set of molecules, a nd 
evolve s stochastically through reactions between these molecules. See 



Wilkinson 



(|2009[ ) and references therein for further background. 

Inference for these models is challenging, as the transition density of the model is 
intractable. However simulation from the models is possible, for example using the 
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Figure 1: Estimated g values from indirect inference (blue) and semi-automatic ABC pilot (black) 
and final (red) analyses, plotted against true g values for 50 g-and-k datasets of sample size n. 



Gillespie Algorithm (IGillespiel . Il977l ). As such they are natural applications for 



ABC methods. Here we foc us on a si mp 
the Lotka-Volterra model of iBoys et al. 



e exa mple of a stoch astic kinetic network, 



(120081 ). While simple, IBovs et al.l tooi) 



show the difficulty of full-likelihood inference. 

This model has a state which consists of two molecules. Denote the state at time t 
by Yt = (y}^\ Y^'^^). There are three types of reaction: the birth of a molecule of 
type 1, the death of a molecule of type 2, and a reaction between the two molecules 
which removes a type 1 molecule and adds a type 2 molecule. (The network is called 
the Lotka-Volterra model, due to its link with predator-prey models, type 1 
molecules being prey and type 2 being predators.) 

The dynamics for the model are Markov, and can be specified in terms of transition 
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over a small time-interval 6t. For positive parameters 6*1, 62 and 6*3, we have 



Pr(Yi+5i = (^i,Z2)|Yi = (2/1,2/2)) 

1 - {9iyi + 92yiy2 + 0^y2)5t 
OiyiSt + o{5t) 
6*22/12/25^ + o{5t) 
e^ysSt + oiSt) 
o{5t) 



o{6t) 



if zi = yi and Z2 = 2/2, 
if zi = 2/1 + 1 and 2:2 = 2/2, 
if Zi = 2/1 - 1 and ^2 = 2/2 + 1, 
if zi = 2/1 and ^2 = 2/2 - 1, 
otherwise. 



We will focus on the case of the network being fully observed at discrete 
time-points, and also on just observing the type 2 molecules initially, together with 
the t ype 1 molecules a t all observation points. All simulations use the parameters 
from bovs etaP hoO^ . with 9i = 0.5, 62 = 0.0025, and 63 = 0.3, and evenly 
sampled data collected at time-intervals of length r. We will perform analysis 
conditional on the known initial state of the system. 

Sequential ABC Analysis 



Wilkinson! (120111 ) consider simulation-based approaches for analysing stochastic 



kinetic networks wh ich are based upon sequential Monte Carlo methods (see 



Doucet et al. 



20o"ol . for an introduction). I n some app l icatio ns, to get these methods 



to work for reasonable computational cost, IWilkinsonI (1201 ll ) suggests using ABC. A 
version of such an algorithm (though based on importance sampling rather than 
MCMC for the sequential update) for the Lotka-Volterra model is given in 
Algorithm m Note there are many approac hes to improy e the c omp utational 



efficie ncy of this algorithm see for example iDoucet et al. 
J2OO1I ) for details. 



(|2000f ) and 



Doucet et al. 



The discussion following Theorem [2] shows that an algorithm like Algorithm H] may 
give inconsistent parameter estimates, even when ignoring the Monte Carlo error. A 
simple remedy to this is to implement noisy ABC within this algorithm. This can 
be done by adding noise to the observed values. The noisy sequential ABC 
algorithm is given in Algorithm O 

To evaluate the relative merits of the two sequential ABC algorithms, we analysed 
100 simulated data sets for the Lotka-Volterra model. For stability of the sequential 
Monte Carlo algorithm we chose K{-) to be the density function of a bivariate 
normal random variable, as a uniform kernel can lead to iterations where all weights 
are 0. We analysed two data scenarios, with r = 0.1, one with full observations, and 
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Sequential ABC for Lotka-Volterra Model 
Input: A set of times, to,. . . ,tn, and data values, yt^ , . . . , y^,^ . 

A number of particles, A^, a kernel, A'(-), and a bandwidth, h. 

Initialise: For i = 1, . . . , A^ sample 6'(*) = {6^^ , 6*^'^ , 6*^'^ ) from the prior distribution for the parameters. 

Iterate: For j = 1, 2, . . . , n 

1. For i — 1, . . . , A^, sample a value for the state at time tj, y^*'', given its value at time tj-i, 
ytj_i, and 0^*^ using the Gillespie algorithm. 

2. For i = 1,. . . ,N, calculate weights 

u;W=i^([y«-y,J//i) 

3. Sample A^ times fro m a Kernel density approximation to weighted sample of 9 values, 



{6iW,mi}£i (see e.g. lLiu and Westl . boOll ). Denote this sample 0^^\ . . . , 9'-'^\ 



Output: A sample of 9 values. 



Algorithm 4: A sequential ABC sampler for the Lotka-Volterra model. 

one where only the number of type 1 molecules is observed. The sequential ABC 
algorithms were implemented with = 5, 000 and h = ^/T. The latter was chosen 
to be a small value for which the sequential algorithms still performed adequately in 
terms of Monte Carlo performance (as measured by variability of the weights after 
each iteration). 

Results are shown in Figure [2j These show both the absolute bias and the root 
mean quadratic loss for estimating each parameter for both ABC algorithms as the 
number of observations analysed varies between 1 and 200. For the full observation 
case, we can see evidence of bias in all three parameters for the standard sequential 
ABC algorithm. Noisy ABC shows evidence of being asymptotically unbiased as the 
number of observations increases. Overall, noisy ABC appears more accurate, except 
perhaps if only a handful of observations are made. When just type 1 molecules are 
observed, the picture is slightly more complex. We observed evidence of bias for the 
standard ABC algorithm for 6i, and for this parameter noisy ABC is more accurate. 
For the other parameters, there appears little difference between the accuracy and 
bias of the two ABC algorithms. The difference between the parameters is likely to 
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Noisy Sequential ABC 

Replace step 2 in Algorithm H] with: 



2. Simulate Xj from A'(x). For i = 1, . . . , iV, calculate weights 



Algorithm 5: Change of Algorithm |4] for noisy sequential ABC. 



be because 6i only affects the dynamics of the observed molecule. 



4.4 Inference for Ricker model 



The Ricker map is an example of an ecological model with complex dynamics. It 
updates a population size Nt over a time step by 



t+i 



rNtC 



-Nt+et 



where the et are independent A^(0, a1) noise terms. IWoodI (120101 ) studies a model in 
which Poisson observations yt are made with means 0iVf. The parameters of interest 
are 9 = (logr, cTf,, 0). The initial state is A'q = 1 and the observations are 
ysi, 2/52, • • • , yioo- This is a complex inference scena rio in w hich there is no obvious 



Wood 



torn argues that 



natural choice of summary statistics. Furthermore 
estimating parameters by maximum likelihood is difficult for this model due to the 
chaotic nature of the system. 



WoodI (120101 ) analyses this model using a Metropolis-Hastings algorithm to explore 
a "synthetic likelihood", Ls{9). A summary statistic function must be specified as in 
ABC, and Ls{9) is then the density of s^j^g under a N{fig, T,g) model where fig and 
Se are an approximation of the mean a nd var i ance of s^^^g for the given parameter 



value obtained through simulation (see 
We simulated 50 datasets with logr 



Wood 



2010 



for more details). 



3.8, = 10 and logcTg values drawn from a 
uniform distribution on [log 0.1,0] (Initial analysis showed that the true log CTg value 
affected the performa nce of t he dif ferent methods). These were analysed under our 
approach and that of IWoodI (l2010l ) (using the code in that paper's supplementary 
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Figure 2: Plots of root square error loss (full-lines) and absolute bias (dashed lines) for the standard 
sequential ABC algorithm (black), and the noisy ABC version (red) as a function of the number of 
observations. Top row is observing the number of both molecules; bottom row is for only observing 
the number of type 1 molecules. Both have t — 0.1. 



material). The distribution just mentioned was used as prior for logcXe, and 
improper uniforms priors were placed on log r > and (p. All parameters were 
assumed independent under the prior. As each parameter had a non-negativity 
constraint, the MCMC algorithms used log transforms of 6 as the state. Each 
semi-automatic ABC analysis used 10^ simulated datasets. 

The summary statistics used by IWoodI (120 lOl ) were the autocovariances to lag 5, the 
coefficients of a cubic regression of the ordered differences At = yt — yt~i on those of 
the observed data, least squares estimates for the model ?/°^^^ = + /32Z/?'^ + e*, 

the mean observation, y, and Y^t^bi^iVt ~ 0) i^^^ number of zero observations). We 
denote this set Eq and use it in the semi-automatic ABC pilot runs. The pilot runs 
used acceptance kernel A = where S is the sample variance matrix of 500 
simulated summary statistics vectors for a representative fixed parameter value. 
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In the regression stage of semi-automatic ABC, training datasets mostly comprised 
of zeros were fitted poorly. Since the observed datasets had at most 31 zeros, any 
datasets with 45 or more zeros was discarded from our training data, and 
automatically rejected in the subsequent ABC runs. Summary statistics were 
constructed for two nested sets of explanatory variables. The smaller set. El, 
included EO and additionally J^tT^i ^iVt = j) for 1 < j < 4, logy, log of the sample 
variance, logX]i=5i ul 2 < j < 6 and autocorrelations up to lag 5. The larger set, 
E2, also added {yt)5i<t<ioo (time ordered observations), {y{t))i<t<50 (magnitude 
ordered observations), (y^ )5i<i<ioo, (2/^j))i<t<50, (log(l + yt))5i<t<ioo, 
(log(l + y{t)))i<t<50, (A?) 52<t<100 

and (A2^)) 

Using E2 instead of El reduced BIC in each linear regression by the order of 
thousands for all but one dataset, suggesting that E2 gives better predictors. Thus 
we used the summary statistics based on E2 within the ABC analysis. The results 
for these ABC analyses show an improvement over the synthetic likelihood for 
estimating logr and 0, and identical performance for estimating Ce- The 
semi-automatic ABC analysis also does better than the comparison ABC one (based 
on summary s tatistics En). No t e that for this application the linear regression 



adjustment of iBeaumont et al.l (120021 ) actually produces worse results then using the 
raw output of the comparison ABC analyses. 

Finally we looked at 95% credible intervals that were constructed from the 
semi-automatic ABC and synthetic likelihood method. The coverage frequencies of 
these intervals were 0.86, 0.70, 0.96 for synthetic likelihood and 0.98, 0.92, 1 for our 
method. Whilst the synthetic likelihood intervals appear to have coverage 
frequencies that are too low for two of the parameters, those from ABC are 
consistent with 0.95 coverage given a sample size of 50 datasets. 



log r CTg 



Synthetic likehhood 0.050 0.032 0.66 

Comparison 0.039 0.038 0.54 

Comparison + regression 0.046 0.041 0.78 

Semi-automatic ABC 0.039 0.032 0.36 

Table 3: Mean quadratic losses for various analyses of 50 simulated Ricker datasets. 
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4.5 Inference for M/G/1 Queue 



Queuing models are an example of stochastic models which are easy to simulate 
from, but often have intractable likelihoods. It has been suggested to analyse such 
models using simulation-based procedures, a nd we will look at a specifi c M/G/1 
queue tha t has been analysed by 
inference (jHeggland and Frigessi . 



3oth ABC (IBlum and Francois 



2OI0I ) and indirect 



20041 ) before. In this model, the service times 



uniformly distributed in the interval [^1,^2] and inter-arrival times exponentially 
distributed with rate ^3. The queue is initially empty and only the inter-departure 
times yi,y2, ■ ■ ■ ,1/50 are observed. 

We analysed 50 simulated datasets from this model. The true parameters were 
drawn from the prior under which (^1, 62 — 61, 6^) are uniformly distributed on 
[0, 10]^ X [0, 1/3]. This choice gives arrival and service times of similar magnitudes, 
avoiding the less interesting situation where all yi values are independent draws 
from a single distribution. 



The analysis of iBlum and Francois! ( 120101 ) used as summary statistics evenly spaced 
quantiles of the interdeparture times, including the minimum and maximum. Our 
semi-automatic ABC pilot analyses replicate this choice, using 20 quantiles. The 
explanatory variables /(y) we used to construct summary statistics were the 
ordered interdeparture times. Adding powers of these values to /(y) produced only 
minor improvements so the results are not reported. The analysis of each dataset 
used 10'' simulated datasets, split in the usual way. 



We also applied the indirect inference approach of iHeggland and Frigessi! ( !2004! ). 
This used auxiliary statistics (|/, min|/j, ^l^'^) where 6^^ is the maximum likelihood 
estimate of 62 under an auxiliary model which has a closed form likelihood; 
independent observations from the steady state of the queue. Numerical calculation 
of 62^^ is expensive so indirect inference used many fewer simulated datasets than 
ABC but had similar runtimes. 

Table !1! shows the results. Semi-automatic ABC outperforms a comparison analysis 
using 20 quantiles as the summary statistics, but once a regression correction is 
applied to the latter the results become very similar. Note that here the 
semi-automatic linear predictors are less accurate when used directly rather than in 
ABC. Indirect inference is more accurate at estimating 62, presumably due to the 
accuracy of 6^^ as an estimate of 62. However it is still substantially less accurate 
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for the other two parameters. One advantage of indirect inference, is that as it 
requires less simulations to estimate the parameters accurately, it can more easily 
accommodate summaries that are expensive to calculate, such as O!^^^. 





Oi 


62 


O3 


Comparison 


1.1 


2.2 


0.0013 


Comparison + regression 


0.020 


1.1 


0.0013 


Semi-automatic ABC 


0.022 


1.0 


0.0013 


Semi-automatic predictors 


0.024 


1.2 


0.0017 


Indirect inference 


0.17 


0.46 


0.0035 



Table 4: Mean quadratic losses for various analyses of 50 M/G/1 datasets. 



4.6 Inference of Tuberculosis Transmission 



Tanaka et al.l ( 120061 ) use ABC to analyse Tuberculosis bacteria genotype data 
sampled in San Francisco over a period from 1991 to 1992. Table [5] shows the data, 
consisting of 473 bacteria samples split into clusters which share the same genotype 
on a particular genetic marker. Thus the data consists of 282 bacteria samples that 
had unique genotypes, 20 pairs of bacteria that had the same genotype, and so on. 



Cluster size 


1 2 


3 


4 


5 8 


10 


15 


23 


30 


Number of clusters 


282 20 


13 


4 


2 1 


1 


1 


1 


1 



Table 5: Tuberculosis bacteria genotype data. 

The proposed model was based on an underlying continuous time Markov process. 
Denote the total number of cases at time t by N{t). The process starts at t = 
with A^(0) = 1. There are three types of event: birth, death (encompassing recovery 
of the host) and mutation. The rate of each event type is the product of N{t) and 
the appropriate parameter: a for birth, 5 for death and 9 for mutation. It was 
assumed that each mutation creates a completely new genotype. Cluster data is a 
simple random sample of 473 cases taken at the first t such that N{t) = 10, 000. 
The model conditions on such a t existing in the underlying process. 

This data contains no information on time, so for A; > 0, parameter values {a, 6, 9) 
and (fca, k5, k9) give the same likelihood. We reparameterise to (a, d, 9) where 
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vrf g, d, 9) oc 7t(9 



m 



Tanaka et al. 



a = a/{a + 6 + 6) and d = 6/{a + 6 + 6). The likelihood under this parameterisation 
depends only on a and d. To reflect prior ignorance of (a, d) we use the prior density 
1 (0 < d < a)I{a + d < 1), where 7r{6) is the marginal prior for 6 used 
(120061 ). The prior restriction d < a avoids the need for simulations 
in which N{t) = 10, 000 is highly unlikely to occur. The other prior restrictions 
follow from positivity constraints on the original parameters. Under this prior and 
parameterisation, the marginal posterior of 6 is equal to its prior and the problem 
reduces to inference on a and d. 



Tanaka et al. 



(120061 ) used two summary statistics for their ABC analysis: and 
H = 1 — '^i{ni/4:73)'^, where g is the number of distinct clusters in the sample, and 
rii is the number of observed samples in the ith genotype cluster. We retain this 
choice for our semi-automatic ABC pilot and the comparison ABC analysis reported 
below. 

As parameters in the pilot output are highly correlated (see Figure [3]), we fitted a 
line to the output by linear regression and made a reparameterisation 
[u vY' = M{a d)'^ where M is a rotation matrix chosen so that on the fitted line u 
is constant. The semi-automatic ABC analysis was continued using {u,v) as the 
parameters of interest. The explanatory variables /(y) comprised: the number of 
clusters of size i for 1 < i < 5, the number of clusters of size above 5, average cluster 
size, H, the size of the three largest clusters, and the squares of all these quantities. 
The semi-automatic ABC analysis used 4 x 10^ simulated datasets in total. 

Figure [3] shows the ABC posteriors, indicating that our methodology places less 
weight on the high d tail of the ABC posterior, reducing the marginal variances 
from (0.0029,0.0088) (comparison) to (0.0017,0.0048) (semi-automatic ABC). 

We further investigated this model through a simulation study. We constructed 50 
data sets by simulating parameters from the ABC posterior, and simulating data 
from the m odel for each of these pairs. We then compared ABC with the summary 
statistics of 



Tanaka et al 



(120061 ) to semi-automatic ABC. Averaged across these 
data sets, semi-automatic ABC reduced square error loss by around 20%-25% for 
the two paramet e rs. H ere we also found that regression co rrection of 

' (l2006h 



Beaumont et al.l (2002) did not change the accuracy of the 



Tanaka et al. 



method. Using the linear predictors obtained within semi-automatic ABC to 
estimate parameters (rather than as summary statistics) gave substantially worse 
estimates, with mean square error for a increasing by a factor of over 3. 
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Comparison 



Semi-automatic ABC 
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Figure 3: ABC output for the Tuberculosis application. Every 1000th state is plotted. The blue 
crosses show the ABC means. 



5 Discussion 



We have argued that ABC can be justified by aiming for cahbration of the ABC 
posterior, together with accuracy of estimates of the parameters. We have 
introduced a new variant of ABC, called noisy ABC, which is calibrated. Standard 
ABC can be justified as an approximation to noisy ABC, but one which can be 
more accurate. Theoretical results suggest that when analysing a single data set, 
using a relatively small number of summary statistics, that standard ABC should be 
preferred. 

However we show that when attempting to combine ABC analyses of multiple data 
sets, noisy ABC is to be preferred. Empirical evidence for this comes from our 
analysis of a stochastic kinetic network, where using standard ABC within a 
sequential ABC algorithm, leads to biased parameter estimates. This result could 
be important for ABC analyses of population genetic data, if inferences are 
combined across multiple genomic regions. We believe the ideal approach would be 
to implement a Rao-Blackwellised version of noisy ABC, where we attempt to 
average out the noise that is added to the summary statistics. If implemented, this 
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would lead to a method which is more accurate than noisy ABC for estimating any 
function of the parameters. However, at present we are unsure how to implement 
such a Rao-Blackwellisation scheme in a computationally efficient manner. 

The main focus of this paper was a semi-automatic approach to implementing ABC. 
The main idea is to use simulation to construct appropriate summary statistics, 
with these summary statistics being estimates of the posterior mean of the 
parameters. This approach is based on theoretical results which show that choosing 
summary statistics as the posterior means produces ABC estimators that are 
optimal in terms of minimising quadratic loss. We have evaluated our method on a 
number of different models, comparing both to ABC as it has been impleme nted in 
the lit erature, to indirect inference, and the synthetic likelihood approach of 



hold ). 



Wood 



The most interesting comparison was with betwee n semi -automatic ABC and ABC 
using the regression correction of lBeaumont et al.l ( 120021 ). In a number of examples. 



these two approaches gave similarly accurate estimates. However, the 
semi-automatic ABC seemed to be more robust, with the regression correction 
actually producing less accurate estimates for the Ricker model, and not improving 
the accuracy for the Tuberculosis model. 

There are alternatives to using linear regression to construct the summary statistics. 
One approach motivated by our theory, is to u se approximate estim ates for each 
parameter. Such an approach has been used in 



Wilson et al. 



(l2009l ). where estimates 



of recombination rate under the incorrect demographic mod el were used as summa ry 
statistics in ABC. It is also the idea behind the approach of iDrovandi et al.l (120 llh . 



where the authors summarise the data through parameter estimates under an 
appro ximating model. A further alternative is to use sliced- inverse regression (jLj, 
199ll ) rather than linear regression to produce the summary statistics. The 
advantage of sliced-inverse regression is that it is able, where appropriate, to 
generate multiple linear combinations of the explanatory variables, such that the 
posterior mean is approximately a function of these linear combinations. Each linear 
combination could then be used as a summary statistic within ABC. 

Finally we note that there has been much recent research into the computational 
algorithms underpinning ABC. As well as the rejection sampling, importance 
sampling and MCMC algorithms we have mention e d, there has been work o n using 



sequential Monte Carlo methods (IBeaumont et al. 



2009 



Sisson et al. 



2007 
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Peters et al. 



20101), and adaptive methods (IBortot et al 



2007 



Del Moral et al 



20091 ) amongst others. The ideas in this paper are focussing a separate issue within 



ABC, and we think the semi-automatic approach for implementing ABC that we 
describe can be utilised within whatever computational method is preferred. 

Acknowledgements: We thank Chris Sherlock for many helpful discussions, and 
the reviewers for their detailed comments. 



Appendix 

Appendix A: Variance of IS- ABC and MCMC-ABC 

Calculating the variance of estimates of posterior means using importance sampling 
are complicated b y the depen dence between the importance sampling weights and 



values of h{6) (see iLiul . 119961 ). A common approach to quantifying the accuracy of 
importance sampling is to use an effective sample size. If the weights are normalised 
to have mean 1, then the e f fectiv e sample size, is divided by the mean of the 



square of the weights. iLiul (119961 ) argue that for most functions h{6), the variance of 
the estimator in (jlj) will be approximately ([5]) but with A^acc replaced by A^g. 

For a given proposal distribution g{6), the effective sample size is 

It can be shown that the optimal proposal distribution, in terms of maximising the 
^eff is 

^opt(^lyobs) °^ ^(^)^^(^|Sobs)^^^- 

In which case 

Var^(p(^|So|3g)) 



where the variance and expectation on the right-hand side are with respect to 7r(6'). 
It is immediate that A'*g > A'acc, with equality only if 'p{Q\?,^^ does not depend on 
Q. The potential gains of importance sampling occur when 'p{Q\&^^ varies greatly. 

Analysis of the Monte Carlo error within the MCMC ABC algorithm is harder. 
However, consider fixing a proposal kernel g{-\-\ which will fix the type of 
transitions attempted. The Monte Carlo error then will be primarily governed by 
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the average acceptance probability. For simplicity assume that K{-) is uniform 
kernel, and that either g{-\-) is chosen to have the prior 7r(^) as its stationary 
distribution, or that the term 7r(^)(yf(6'j_i|6')/{7r(^^i_i5f(6'|6',j_i)} ^ 1 and can be 
ignored. The average acceptance probability at stationarity is 

1 1 T„o(»|Sobs)9(»'l»)P(»'Nobs)d»d#' = ( j ir(»)p(»|Sobs)d»') x 

The integral comes from averaging over the current and proposed values for the 
MCMC algorithm, with the average acceptance probability for a given proposed 
value 6' being p{6'\Sq^^). The right-hand side comes from using ([1]). The first term 
on the right-hand side is the average acceptance probability of the rejection 
algorithm. The second term is 1 if we use an independence sampler g{6'\6), and will 
be >> 1 if the ABC posterior is peaked relative to the prior, and if the transition 
kernel proposes localised moves. 



Appendix B: Proof of Lemma [T] 

Write 

pie\s^^^)ni9)de = J J Kis-s^^^)/h)nis\9)nie)deds 



j /i'^K(x)7r(sQ]3g + /ix)dx. 



The first equality comes from the definition of p{6\Sq^^), the second by integrating 
out 6 and making a change of variable x = (s — s^j^g)//?,. So 

h-'^ J p{9\s^^^)iT{9)d9 - 7i{s^l^^) < y"fs:(x)|7r(sQbg + /ix) -7r(Sobs)|dx 
This bound is shown to be o(l) under either condition. 

Consider first condition (i). Define c to be the maximum value of |x| such that 
J^dxl) > 0. For any e > 0, by continuity of -7r(s) at s^j^g we have that there exists a 
6 > such that |x| < 6 implies ^(sQi^g + x) — 7r(sQ]^g)| < e. Define he = S/c. Then 
for h < hf we have 



//f(x)|.(s„b, + ftx)-.(s„b,)|dx<e. 
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This inequality follows as h\x.\ < 6 for all x where K[x.) > 0. 

Now consider condition (ii). By differentiable continuity, Taylor's theorem gives 



TT 



The remainder factor rj(x) is \d7r{z)/dsi\ for some z(x), so by assumption, 
< R, a. finite bound. Thus 

y _ft'(x)|7r(sQ]^g + /ix) — 7r(sQ]gg)|dx < hR'^^ J \^i\K{x)dx. 



Appendix C: Proof of Theorem [4] 

Rearrangement of the loss function gives 

EiL{9, e- A)|yobs) - trace(AS) = E((e~ - efA{e - e)\y^^^), 

where 6 = E^Oly^^^), the mean under the true posterior, which equals 'S'(yQ]^g) for 
the summary statistics under discussion. Define "^(sQi^g) — ^(^obs) — ^obs '^^^ make 
the change of variables x = (6' — s^^^)/h. Now 

E{L{9, 9; A)|yobs) " trace(AS) - j x^Axfs:(x)dx 

= -2h j i^^A5{e + /ix)/s:(x)dx + j ^{d + hxfAS{9 + /ix)/s:(x)dx. 

It is required that the modulus of the right hand side is o{h'^). Let R be the support 
of K. Since this is finite, it suffices to show that 

max5(^ + hx.) = o{h). 

X6-R 

To find an expression for (5, observe that ([I]) holds for noisy ABC. Taking its 
expectation and making the change of variable s = S'(y) gives 

. ^ jeix{e)ix{s\e)K{{s-t}/h)deds 

^ ' /7r(s)A'({s-t}//i)ds 

Note 

97r{9)7r{s\9)d9 = 7r(s)E(^|s) = S7r(s 
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where the second equahty is due to our choice of S. Thus 

^ /(s-tMs)K({s-t}//^)ds 
^ ' /7r(s)Ji({s-t}//i)ds ■ 

Make the change of variables y = (s — t)/h and consider the case t = 6 + /ix, then 

hJy7r{e + h[^ + y])K{y)dy 



5{e + hx.) 



J n{9 + h[^ + y])K{y)dy 



By the argument in Appendix B, continuity of 7r(s) at s = 6' gives denominator 
7t{6) + 0(1). Consider the ith component of the integral in the numerator 



y,7r{e + /i[x + y])K(y)dy - j y,7i{e)K{y)dy 



< max l-u,- 
- y&R 



V(^~ + /,[x + y])-7r(^~))Aly)dy 



This integral is o{h) by the continuity argument just mentioned. Noting that 
/ yiK{y)dy = by assumption, we have 



I- 



y7r(^ + /i[x + y])J^(y)dy = 0(1). 
Combining these results gives the required bound for 5. 
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